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COMPARING THE NEATM WITH A ROTATING, CRATERED THERMOPHYSICAL ASTEROID MODEL 
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ABSTRACT 

A cratered asteroid acts somewhat like a retroflector, sending light and infrared radiation back toward the 
Sun, while thermal inertia in a rotating asteroid causes the infrared radiation to peak over the "afternoon" part. 
In this paper a rotating, cratered asteroid model is described, and used to generate infrared fluxes which are 
then interpreted using the Near Earth Asteroid Thermal Model (NEATM). Even though the rotating, cratered 
model depends on three parameters not available to the NEATM (the dimensionless thermal inertia parameter 
and pole orientation), the NEATM gives diameter estimates that are accurate to 10% RMS for phase angles less 
than 60°. For larger phase angles, such as back-lit asteroids, the infrared flux depends more strongly on these 
unknown parameters, so the diameter errors are larger. These results are still true for the non-spherical shapes 
typical of small Near Earth objects. 
Subject headings: asteroids, size, infrared 



1. INTRODUCTION 

Rotation causes a diurnal oscillation in the illuminating flux 
on a surface element of an asteroid. During the day, heat is 
conducted into the surface, while during the night this heat is 
radiated. The combination of this phase delayed conducted 
heat and the direct heat from the Sun leads to a tempera- 
ture maximum during the afternoon on an aste roid, just as i t 
does on the E arth. Rotating m odels of planets dWrightll 19761) 
and asteroids dPeterson II 1976b which incorporate the effects 
of thermal inertia have been in use for decades. But these 
models do not show the peaking near zero phase angle seen in 
real asteroids. This lack is addressed in the stand ard thermal 
model (STM) for asteroids dLebofskv et al.ll 19861) by evaluat- 
ing the flux at zero phase angle and then applying a linear 0.01 
mag/degree phase correction. This beaming of the infrared ra- 
diation toward the Sun reduces the total reradiation, so in or- 
der to conserve energy the subsolar temperature is computed 
by replacing the emissivity of the surface e by e-q, where the 
beaming correction 77 = 0.756. This approximately conserves 
energy, but the STM is really an empirical fitting function 
rather than a physical model. Infrared observations of Near 
Earth Objects (NEOs) analyzed using the STM yielded inac- 
curate diameters for small, rap idly rotating NEOs seen at large 
phase angles, so lHarrisI d 19981) developed the Near Earth As- 
teroid Thermal Model (NEATM). In the NEATM the beam- 
ing correction 77 is an adjustable parameter, and the infrared 
flux is evaluated by integrating the emission over the asteroid 
surface seen from the actual position of the observer. The ob- 
served color temperature is matched by adjusting rj, and then 
the NEATM specifies the average surface brightness of the 
object, so the observed flux implies a diameter. 

Neither the STM nor the NE ATM has a phy sical explana- 
tion for the beaming effect, but [Hansen] d 19771) provides one 
using a cratered asteroid model. If an asteroid is covered with 
craters, the peak temperature on the surface at the sub-solar 
point will be higher. Furthermore, the flux at zero phase an- 
gle from areas near the limb will be much higher than in ei- 
ther STM or the NEATM, since while the surface is a mixture 
of illuminated and shadowed areas, the illuminated areas are 
closer to facing the Sun, and an observer at zero phase angle 
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sees only the lighted areas. The lack of visible shadows at 
zero phase angle cause a peak in the emission at opposition. 
In order to get enough shadowing, it is necessary to have a 
substantial amount of concavity in the surface, which could 
be due to a porous gra nular surface, or to the craters consid- 
ered bv lHansenl d 1977b . 

Cratering and thermal inertia have been comb ined into a 
single thermophysical model by iLagerrosI (1 19961) . The ther- 
mophysical model used in this paper has been develo ped in- 
dependently, but is very similar to the lLagerrosI 0996) model 
with the crater covered fraction set to 100%, and with depth 



to width ratio given by ( 1 
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2. ROTATING CRATERED MODEL 




FIG. 1. — Physical analog of the cratered asteroid model, seen at ~ 60° 
phase angle. In the actual computational model the flat areas between the 
craters are eliminated. 
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FIG. 2. — Heat flow into the surface (G) for a triangle temperature forc- 
ing function T. Computed with 32 samples per period. Two full cycles are 
shown. 




FIG. 3. — Temperature distribution in the late afternoon for a crater at 
4/ circ latitude when the sub-solar latitude is 30°. The thermal inertia pa- 
rameter is = 1. The Sun is setting in the Northwest, so the Southeast rim 
of the crater is still sunlit. 
-11 

' ' I 



10" 



10 



16 



"T~ 



"T~ 



T 



T 



T 




20 40 60 80 100 
Phase angle [degrees] 



120 



FIG. 4. — Infrared phase curves for an asteroid 1.4 AU from the Sun in the 
passbands of the IRAC on the Spitzer Space Telescope, 3.6, 4.5, 5.6 & 8 fim 
from bottom to top. The solid curves show a cratered model with maximum 
slope 9 max = 45°, while the dotted curves show an uncratered model. The 
dashed line shows the 0.01 mag/degree phase curve of the STM. 



In the rotating cratered asteroid model, heat is conducted 
vertically into and out of the surface, but not horizontally. The 
facets in a given crater can see each other and see the Sun, but 
there is no radiative transfer of heat from one crater to another. 
Thus the problem of finding the temperature distribution on 
the asteroid breaks up into many small problems of finding 
the temperature vs. time of each facet. 

The equation describing heat conduction in the surface 
layer is 

where k is the thermal conductivity, p is the density, and C is 
the specific heat per unit mass. The heat flow into the surface 
is Letting Q = z/n, then the heat flow is ^ and the 

heat conduction equation is 

upC 



dT 
~dt 



d 2 T 



(2) 



which depends only on the thermal inertia T = y 'npC. 
The units of T are [E / \t(T '/ ' L)L 2 )M / X 3 E / '{TM)] 1 ' 2 = 
E/{TL 2 Vt). In IWrightl (fl976b the value of T = 
0.006 cal/cm 2 /K/sec 1//2 was use d for Mars. In more modern 
units this is 251 J/m 2 /K/sec 1/2 . lHarris I (120061) gives values 
for r of 10-20 for main belt asteroids, 50 for the Moon, 150 
& 350 for the NEOs Eros & Itokawa, and 2500 J/m 2 /K/sec 1 / 2 
for bare rock. 

Putting in the solar heating and thermal radiation, I get 
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with 6 = fit where il is the rotational frequency of the as- 
teroid, R is the distance of the asteroid from the Sun, Lq is 
the solar luminosity, A is the albedo, e is the emissivity in the 
thermal infrared, and a is the Stephan-Boltzmann constant. 
This assumes the special case of an equatorial region with the 
sub-solar latitude (/?q) of zero. In general the max(0, cos#) 
would be replaced by either the cosine of the angle between 
the Sun and the surface normal, or zero if the facet is shad- 
owed on the Sun is below the horizon. 
The unit of temperature in the code is 
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which is the equilibrium temperature of a surface oriented to- 
ward the Sun. Let y = T/T , giving 



and 
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Redefine the depth variable again using w = eaT^ = z/H 
where H is the distance such that a temperature gradient of 
T /H sets up a conductive flux equal to the radiative flux 
eerTo . These equations are now 
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FIG. 5. — Temperature distribution as a function of latitude and local time for a rotating cratered model with thermal inertia parameter © = 0.1 and maximum 
crater slope 9 ma x = 45°. The sub-solar latitude is 30°. 
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FIG. 6. — Temperature distribution as a function of latitude and local time for a rotating cratered model with thermal inertia parameter 
crater slope 9 m ax = 45°. The sub-solar latitude is 30°. The circle just above the equator at local time 17:37 is blown up in Figure]!] 
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is a dimensionless measure of the importance of thermal in- 
ertia on the temp eratures. It is a bout = 1.5 for Mars using 
the parameters of lWrightl (119761) . 
Now let y = J2 n Vn exp(wi0 — k n w). The solutions to the 



must have inQ 2 = k\ so k n = 0-^/n/2(l + i) for n > 0. For 

negative n, k n = 9i/|n|/2(l — i) must be taken to guarantee 
the solution is damped toward z = oo. 

The physical length scale given by H/Q = T/(CpVQ). 
Since Cp w 10 6 J/m 3 /K and O w 1CT 4 rad/sec, the length 
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FIG. 7. — Temperature distribution as a function of latitude and local time for a rotating cratered model with thermal inertia parameter = 10 and maximum 
crater slope 8 m ax = 45°. The sub-solar latitude is 30°. 
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FIG. 8. — Definition of the clock angle. marks the sub-solar point, and © 
marks the sub-observer point. The angle between the pole and the sub-solar 
point is 90 — Pq . 

scale is smaller than 10 cm. Thus the assumption of no hori- 
zontal heat conduction is reasonable for craters > 1 m. 
The heat flow into the surface is 

<9y = _y- 
dw ^ 




y n k n exp(in6>). 



(11) 



Now take y(6, 0) as a triangle wave centered at 6 = with 
y = 1, and linearly sloping down to y = at 6 = ±2ir/N. 
This gives 



y n = Tr 



2tt/N 
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FIG. 9.— Flux in the WISE 12 and 23 /im bands for a 1 km asteroid 1 AU 
from the Earth and 1 .4 AU from the Sun, observed at 50° phase angle. The 
sub-solar latitude is 30° and the thermal inertia parameter is = 1. The 
clock angle rotates the line of sight around the Sun line while keeping a con- 
stant phase angle. For 0° clock angle the observer is East of the Sun, giving 
an afternoon view. For 90° clock angle the observer is over 80° latitude at 
noon. These fluxes lead to a beaming parameter r\ and calculated diameter D 
using the NEATM. 



The heat flow into the surface is QG(6) with 
7V(1 - cos(2irn/N)) 
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G{6) = ]T 
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The value of this G{6) at 9 = 2nm/N defines the vector G 
used in code. The average of this over bins of width 2-k/N in 
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Figure |2] shows the function G(m) for TV = 32. 

The mutual irradiation of the crater facets introduces a cou- 
pling matrix. The contribution of facet j to facet i goes like 
n~ 1 AQjR 2 cos(8j) cos(0)/cff - where 0j = 8j are the angle 
of emission and incidence, Aftj is the solid angle of the facet 

on the spherical cap crater, and c£y = 2Rcos(dj) is the dis- Averaging this equation gives 
tance between the facets. Thus the coupling is just AO/4.7T. 
The total light falling on a facet is given by 



Di was computed as fj,if v times F@, where /„ is the fraction 
of the facet that is visible from the Sun, and fx is the cosine of 
the angle between the surface normal and the Sun. The visible 
fraction is computed using a fi ner pixelization of the sphere, 
HEALpix (iGorskiet al.ll2005l) with 49,152 pixels. Figure 
shows the facet structure with a crater. 

In the thermal infrared, the mutual visibility of the facets 
couples their temperatures together. The total infrared flux 
falling on a facet is 



in 



+ (1 - e)f c (H) 
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where Di is the direct solar flux, and Si is the total flux on a 
facet. Averaging this equation gives 
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where f c is Sl crater /47r, the fraction of the sphere included 
in the spherical cap craters. For the 8 max = 45° used here, 
f c = (1 — cos8 max )/2 = 0.15. Therefore 

^ = A + ^44 (18) 



1 - Af c 



For the calculations reported here, the craters were divided 
into 127 facets, consisting of a central circle surrounded by 
rings of 6, 12, ... , 36 square facets. Since the resulting facet 
size of 28 max /i3 f=a 7° was fairly coarse, the direct insolation 
10 r 
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FIG. 10. — As in Figure [9] but at 90° phase angle. For 0° clock angle the 
observer is over the equator at 6 PM. For 90° clock angle the observer is over 
60° latitude at midnight. 



H l= <H>^ 
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Remembering that the unit of flux is both (1 — A)Fq and 
eaT^, and that only a fraction e of the incident heat H is ab- 
sorbed, we get a power balance equation 



Vij = (fvfJ-h + y 
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where j/y is the temperature of the i facet at the j time, 
and there are n facets and N times. This is a set of n x N — 
127 x 32 = 4064 coupled non-linear equations in 4064 vari- 
ables. Fortunately the facet to facet coupling is weak and can 
be handled by iteration, so the task of solving for the temper- 
atures is not too onerous. 

There is little point in using extremely fine subdivisions of 
the asteroid's surfa ce. The bright ness temperatures of Mars 
were calculated by I Wright! (1 19761) using only 6 latitudes. In 
the calculations reported here, 16 latitudes were used, and 
the rotation period was divided in 32 time steps. This gives 
65,024 temperatures to be found. These temperature pattern 
are plotted in Figures HD & □ for 6 values of 0.1, 1 & 10 
with sub-solar latitude of 30°. 

2.1. Observed Flux Calculation 

Given the temperature distribution, the observed flux is 
found by integrating over the surface of the asteroid. For a 
given frequency v, the quantity x = hv/kT is found. Then 
the observed infrared flux is found using 



2div z 



^max(/i/(/3,j),0) 
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Ei(/«M) 



<jr> 



(22) 



where dE(/3, j) is the surface area of the asteroid in the bin at 
latitude /3 and longitude given by j, D is the distance to the 
observer, /x/ is the cosine of angle between the normal to the 
surface and the line of sight, and (fvl^ijp is the fraction of the 
i th facet visible by the observer times the cosine of the angle 
between the facet normal and the line of sight. The observed 
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FIG. 11. — The beaming parameter r\ and the fitted diameter as a function of the phase and clock angles for the asteroid in Figure|9] 

This can be multiplied by F u (&)/ Ftoi(&) to give the reflected 
optical spectrum. 

In this paper the latitudes are uniformly spaced in sin (3, so 
dH is a constant for a sphere. The assumption that the as- 
teroid is spherical only enters into dY,, so any convex shape 
for an asteroid can be accommodated merely by changing the 
weights dYj going into the flux sum. It is important to remem- 
ber that the angles (3 and A refer to the surface normal vector, 
not the vector from the center to a surface element. Thus at 
the end of the major axis of a 2 : 1 : 1 ellipsoid, ri£ is a mini- 
mum, 16 times smaller than the value at the end of the minor 
or intermediate axes. 

3. NE ATM FITTING 

The rotating cratered asteroid model described model has 
been used to predict infrared fluxes for NEOs, and then these 
fluxes have been used in the NEATM to find the beaming pa- 
rameter ?y and the diameter d. For a given observation, the 
distance to the Sun and the phase angle are known, but there 
are still two angles and the thermal inertia parameter 9 that 
need to be specified. Figure[8] shows the definition of the two 
angles, which are the sub-solar latitude and the clock angle. 
Examples are shown in Figures l9l and [TOl for a 1 km diame- 
ter asteroid 1.4 AU from the Sun, with albedo A = 0.1 and 
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FIG. 12. — Histogram showing the frequency of various diameter discrep- 
ancies when using the NEATM to analyze fluxes from the rotating cratered 
model, assuming spherical asteroids. 



bolometric optical flux is 



F opt = ^2 max(/z/(/3, j), 0) 
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(23) 
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FIG. 13. — RMS diameter errors for a sample of 576 NEOs that could be observed by the Spitzer Space Telescope in Cycle 4. While the distance from the the 
Sun and phase angle are held fixed at the actual values, fluxes were simulated for Monte Carlo asteroids with random thermal inertias uniform in pole position 
and log © in the range 0.1 < < 10. Fluxes at the WISE 12 & 23 /im bands gave NEATM RMS diameter errors shown by the lower dots. The IRAC 5.8 and 
8 (im fluxes give the upper dots. 
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FIG. 14. — The egg shaped example. The figure is rotationally symmetric 
about the horizontal axis, while the object rotates about a vertical axis. 
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FIG. 16. — The variations of the projected area in blue, the square of the 
NEATM derived diameter in red, and the NEATM r] parameter in green. 
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FIG. 15. — Lightcurves for the egg shaped asteroid seen at 45° phase 
angle with the rotation axis normal to the Earth-Sun-asteroid plane. The solid 
curves show viewing the afternoon side, while the dashed curves show the 
morning side. Black is optical, blue is 3.6 fim, green is 5.8 fim, and red is 12 
fjm. 

emissivity e = 0.9. At a 50° phase angle, the errors in the 
NEATM calculated diameters are small, but for 90° the maxi- 
mum error increases by a factor of 4. These calculations have 
been done using the planned 12 & 23 fim bands of the Wide- 




FlG. 17. — Histogram showing the frequency of various diameter discrep- 
ancies when using the NEATM to analyze fluxes from the rotating cratered 
egg-shaped model. The solid histogram shows the distribution of diameter 
errors derived from single observations at a random rotational phase. The 
dashed histogram shows the distribution of diameter errors derived from the 
average over a rotation. The dotted curve shows a Gaussian with a 10% stan- 
dard deviation. 
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FIG. 18. — RMS diameter errors for a sample of 576 NEOs that could be observed by the Spitzer Space Telescope in Cycle 4. While the distance from the the 
Sun and phase angle are held fixed at the actual values, fluxes were simulated for Monte Carlo asteroids with random thermal inertias uniform in pole position 
and log in the range 0.1 < 010. Fluxes at the WISE 12 & 23 fim bands were used in the NEATM. RMS diameter errors when a single observation is made 
of each asteroid give the upper dots, while errors after averaging the NEATM diameter give the lower dots. 



field Infrared Survey Explorer (WISE, 
which is scheduled for launch in 2009 



dMainzer et alj|2006l) ) 



4. NEATM ACCURACY 

The NEATM can reproduce the model diameters for quite 
well for phase angles a < 60°. Figure [Til shows that the er- 
rors are small for any clock angle as long as the phase angle 
is smaller than 60°. This conclusion does not depend strongly 
on either the sub-solar latitude or the thermal inertia parame- 
ter O. To show this for a representative sample of real obser- 
vations, Monte Carlo simulations of NEO observations from 
a Spitzer Space Telescope proposal (SNEAS, PI Eisenhardt). 
576 NEOs were found to be observable by Spitzer during Cy- 
cle 4, with good signal-to-noise ratio in the IRAC 5.8 & 8 
/mi bands. Only the distance to the Sun and the phase an- 
gle were taken from this observation table. Then fluxes were 
computed for a random distribution of pole positions and ther- 
mal inertias. The pole positions were chosen uniformly in Air 
steradians, and the thermal inertias were chosen uniformly in 
the logarithm in the range 0.1 < 9 < 10. To choose a ran- 
dom pole position one picks sin (3q uniform in [0 — 1] and the 
clock angle uniform in [0— 2n]. For each of the 576 objects 30 
different choices of pole position and O were analyzed. The 
fluxes were then analyzed using the NEATM to derive r\ and a 
diameter. The RMS diameter errors, binned by phase angles, 
are shown in Figure [T_3] Both the WISE 12 & 23 /im bands, 
and the Spitzer IRAC 5.8 & 8 /im bands give data that work 
well with the NEATM. 

When the NEATM breaks down at large phase angles, the 
estimated diameter is usually too large. Figure [121 shows the 
distribution of the errors for the Monte Carlo observations in 
FigureQj] The distribution is clearly positively skewed. 

Figure Qj] shows that the NEATM works reasonably well 
for moderate phase angles when compared to a more com- 
plete thermophysical model. But to show that the NEATM 



works in the real world one needs comparisons to real ob- 
jects with size de t ermine d by radar or spacecraft imaging. 
Harris & Lagerrosl d2002l) find that NEATM diameter errors 
average less than 10% for phase angles less than 60°, but the 
number of objects in the comparison was quite small. 

5. ASPHERICITY 

All of the previous caculations have assumed spherical ob- 
jects, but small NEOs usually have quite aspherical shapes. 
An example of a non-spherical shape is the egg-shaped object 
seen in Figure [14] This is half of a 2.5:1:1 ellipsoid joined 
to half of a 1.25:1:1 ellipsoid. If this object is 1 AU from 
the Earth, and 1.414 AU from the Sun, with phase angle of 
45° and a sub-solar latitude of 0°, thermal inertia parame- 
ter = 1, albedo of 0.1 and emissivity of 0.9, one gets 
the lightcurves seen in Figure [15] for clock angles of 0° and 
180°. The flux normalization applies to a 1 km short axis and 
a 1.875 km long axis. The optical lightcurve amplitude is 1.4 
mag peak to peak while the 23 /im amplitude is 0.84 to 0.93 
mag. These amplitudes are larger than the 1 .875: 1 variation in 
the projected area for this shape. The NEATM applied to the 
12 and 23 fim data gives the rj and D values plotted in Figure 
[161 along with the projected area as a function of rotational 
phase. The diameter from the NEATM tracks the projected 
area fairly well, and there is a definite variation of the beam- 
ing parameter rj with rotational phase. 

The calculations of the accuracy of the NEATM for the 
Spitzer sample of 576 NEOs have been repeated for the egg- 
shaped asteroid shown in Figure [14] A non-spherical shape 
introduces another parameter, the rotational phase, that must 
be either treated as a random variable or integrated over. Fig- 
ure[l7]compares NEATM diameters to the true "diameter" of 
the egg. There are many ways to define the true diameter: 
the diameter of the sphere with the same volume of the egg is 
1.233 km, while the sphere with the same surface area as the 
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egg has a diameter of 1 .275 km. Since the NEATM is trying to 
estimate the projected area of an object it seems reasonable to 
use the equivalent area diameter of 1 .275 km as the reference. 

Even a single sample at a random rotational phase gives 
a reasonable diameter estimate: the ratio of the NEATM to 
equivalent area sphere has a median and interquartile range 
of 1.042^0 096- Since the interquartile range only contains 
50% of the sample, these errors should be considered to be 
standard errors instead of standard deviations. One can do 
better using the NEATM diameter averaged over the rotation 
period. This gives a median ratio and interquartile range of 
1.01 9^0 072- The improvement coming from better sampling 
of the rotational phase is fairly small, even though the egg 
shape considered here gives high peak to peak amplitudes up 
to 1.4 mag in the optical and short infrared wavelengths. The 
small amplitude in diameters is due to several factors: the 
lightcurve amplitude is lower at the WISE wavelengths of 12 
and 23 fim, the rms of a sine wave is 2.8 times less than the 
peak to peak, the amplitude goes down to zero for pole-on 
orientations, and the diameter varies like the square root of the 
flux. The median and interquartile range of the rms variation 
of the NEATM diameter with rotational phase is only T.GjJ'g 
percent, while the median peak to peak 12 fim flux amplitude 
is 0.647 mag. 

It is better to evaluate the diameter at each rotational phase 
using the NEATM and average these diameters than to av- 
erage the fluxes and use these averages in the NEATM. The 
RMS errors for both single sample and rotationally averaged 
NEATM diameters for the sample of NEOs observable by 
Spitzer are shown in Figure[l8] For each asteroid, 30 different 
pole positions and inertia parameters 9 have been simulated, 



using the egg-shaped model discussed above. 

6. CONCLUSION 

The errors in diameters computed from the NEATM for as- 
teroid observations at phase angles less than 60° are less than 
10% RMS, even for the non-spherical shapes typical of NEOs. 
For WISE, observing at 90° elongation, any object with a dis- 
tance larger than 0.6 AU will have a phase a < 60°. The 
Spitzer Space Telescope can observe at elongations between 
85° and 120°, so only very close passes involve a > 60°. 
The error evaluated in this paper only includes the errors due 
to not knowing the thermal inertia and pole orientation of an 
asteroid. There will be additional errors due to uncertainties 
in the true emissivity of the asteroid surface, but these errors 
should be small, since the emissivities in the thermal infrared 
are quite close to the maximum possible value of 1 .0. Diam- 
eters from the NEATM do not depend on the assumed albedo 
so there is no additional error from albedos. WISE will obtain 
« 10 observations of each asteroid spread over 30 hours, and 
will thus get good sampling of asteroid lightcurves, which re- 
duces the errors associated with non-spherical shapes. WISE 
will be sensitive enough to measure hundreds of thousands of 
asteroids, and fitting the WISE 12 & 23 fim fluxes using the 
NEATM will provide reasonably good diameter limits for a 
large sample of asteroids. 
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